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Abstract 

This paper presents an approach for obtaining accurate interaction ener- 
gies at the DFT level for systems where dispersion interactions are important. 
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This approach combines Becke and Johnson's [J. Chem. Phys. 127, 154108 
(2007)] method for the evaluation of dispersion energy corrections and a Hirsh- 
feld method for partitioning of molecular polar izability tensors into atomic con- 
tributions. Due to the availability of atomic polarizability tensors, the method 
is extended to incorporate anisotropic contributions, which prove to be im- 
portant for complexes of lower symmetry. The method is validated for a set 
of eighteen complexes, for which interaction energies were obtained with the 
B3LYP, PBE and TPSS functionals combined with the aug-cc-pVTZ basis set 
and compared with the values obtained at CCSD(T) level extrapolated to a 
complete basis set limit. It is shown that very good quality interaction energies 
can be obtained by the proposed method for each of the examined functionals, 
the overall performance of the TPSS functional being the best, which with a 
slope of 1.00 in the linear regression equation and a constant term of only 0.1 
kcal/mol allows to obtain accurate interaction energies without any need of a 
damping function for complexes close to their exact equilibrium geometry. 

Keywords: Interaction Energies, Dispersion, Hirshfeld, DFT, Basis Set 
Extrapolation. 
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2 Introduction 

Since Density Functional Theory has been introduced into the world of computational 
chemistry, countless studies on "large" systems have been performed that would oth- 
erwise have been impossible due to the size of the systems examined.^ A large portion 
of those studies involve biologically active molecules, their structure, reactivity, cat- 
alytic and binding propertiesP One of the most fundamental aspects examined in 
these studies are interaction energies. However, the use of DFT becomes problematic 
when energetics and related properties are examined for systems where dispersion 
interactions are important P Accurate description of interaction energies demands 
the use of levels of theory that include electron correlation, and although MP2 has 
started to become applicable to systems of relevant size in recent yearsP® the more 
accurate methods such as CCSD(T) are still far from reaching that stage. Therefore, 
the adjustment of DFT methods for a correct descriptions of dispersion interaction is 
nowdays a topic of an active research. Numerous examples can be found in the recent 
special issue of Phys. Chem. Chem. Phys. dedicated to stacking interactionsP 

In several studies adjustments were made to existing density functionals to im- 
prove their performance for non-covalent interactions. Xu et alP designed an X3LYP 
extended functional, based on the well-known B3LYP functional, that improves the 
accuracy for Van der Waals complexes. Zhao and Truhlar^^ developed functionals 
based on simultaneously optimized exchange and correlation functionals. Rothlis- 
berger et al. developed a dispersion-corrected DFT, where they augment the B3LYP 
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functional with dispersion corrected atom-centered potentials (DCACPs).^^ Also 
several studies by Hirao and Lundqvist have been performed for developing special 
correlation functionals which take long-range dispersion interactions into account P2HH1 
The inclusion of empirical dispersion coefficients,^^ as advocated by Grimme, Hobza 
and Head-Gordon, has booked some success as a cost-efficient method for examina- 
tion of stacking phenomena in larger systems, since such applications are today not 
yet possible with the more elaborated methods mentioned above. 

A different approach consists of calculating dispersion energies from dispersion 
coefficients. For instance, Van Gisbergen et al. derived van der Waals dispersion co- 
efficients from frequency dependent polarizabilities using time dependent DFTpiEU 
On the other hand, Becke and JohnsorpDES] developed an approximation for the cal- 
culation of dispersion coefficients from exchange-hole dipole moments that allows to 
obtain dispersion energies at the DFT level in an easy and efficient fashion. However, 
the method of Becke and Johnson has two major drawbacks. First of all, Becke and 
Johnson use inappropriate values for the polarizabilities of atoms in the molecules, 
thus undermining the theoretical foundation of the method significantly. In their 
first publication they used empirical values for polarizabilities, while in their later 
works they approximated atomic polarizabilities by scaling free-atom polarizabilities 
by the Hirshfeld effective volume of the atom in the molecule. Such a rough approach 
for obtaining atom-in-molecule polarizabilities ignores many relevant effects, such as 
electron density reorganization due to the applied electric field. Second, Becke and 
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Johnson make use of a damping function that strongly reduces the values of dispersion 
energy even at equilibrium geometries by an approximate factor of 2. 

In our previous worlP^we showed that the use of intrinsic polarizabilities, obtained 
from the Hirshfeld methodj^SlGSE] improves the dispersion energies obtained from the 
dispersion coefficients of Becke and Johnson significantly. As a result, not only the 
dubious character of the atomic polarizabilities used in the method is eliminated, but 
also realistic dispersion energies are obtained at equilibrium geometry without the 
need for a damping function. 

In this work we develop our all-Hirshfeld approach further in three separate ways. 
First of all, we extend the methodology for reproducing high-level interaction ener- 
gies at the DFT level, instead of comparing pure dispersion energies. This method 
is applied using three functionals of a different nature, namely the hybrid functional 
B3LYP,E2the GGA functional PBEPEIand the meta-GGA functional TPSS, 45 none 
of them including non-local correlation, in order to test its universal character. Sec- 
ond, in order to improve the accuracy for complexes of reduced symmetry, we intro- 
duce anisotropy for the derivation of the dispersion coefficients. Finally, we introduce 
the iterative Hirshfeld method (Hirshfeld-I)0S! into the calculation of the dispersion 
coefficients. The Hirshfeld-I method, recently developed by Bultinck et al.,^ brings 
several fundamental improvements to the classic Hirshfeld weight function and the 
resulting partitioned properties such as charges, dipole moments and polarizabilities. 
By its iterative nature the method eliminates the somewhat arbitrary nature of the 
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weight function and can also be applied to charged systems. It has recently been 
shown by some of the authors that not only the charges obtained with Hirshfeld-I 
are more in line with the oxidation state of the atoms in the molecules, but also 
the intrinsic polarizabilities are more adequate.^ For the validation of our modified 
approach a set of 18 complexes is examined. To ensure the quality of the high level 
geometries and interaction energies, all the complexes were optimized using the same 
methodology, namely CCSD(T)/aug-cc-pVTZ and a frozen monomer approach. 



3 Method 

The dispersion corrections used in this work are based on the model developed by 
Becke and Johnson EHUD wherein explicit expressions for the dispersion coefficients Cq, 
Cs and Cm were derived from the instantaneous dipole moment created by an electron 
and its corresponding Fermi hole. The dispersion energy between two nonoverlapping 
systems A and B at a distance R from each other is given by 

disp ~ \R* + R* + Rio)' [ ' 

According to Becke and Johnson's model the coefficients in eq. ([T]) can be obtained 

from the polarizabilities a of the systems and the expectation values of the square of 

their dipole moments (Mf), quadrupole moments (M|) and octopole moments (Mf) 

a A a B (Mp A {Ml) B 
° 6 a A (Mf) B + a B (Mf} A > {) 

r = 3 olaolb {{Ml) A (Mj) B + (Mj) A (Mp B ) 
8 2 a A {Ml) B + a B {Ml) A ' W 
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r Q ^AOB {{Ml) a (Mj) B + (Mj) A (Mf) B ) 

10 a A (Mt) B + a B (Mt) A 

, 21 olaclb ((Mp A (M*) B + (M|) A (M*) B ) 

5 a A {Ml) B + a B (Mt) A ' 1 ] 

The expectation values of the square of a multipole Mi are approximated by 

( M i) = J2 / Mr)[r'-(r-cW] 2 d 3 r, (5) 

with (j representing the spin of the electron. 

It has also been shown by Becke and Johnson that when the systems A and B 
contain more than one atom, the dispersion energies obtained from coefficients in 
equations fl2JH]) can be decomposed into pair-wise atom-atom interactions between 
the atoms in the two systems. For example, the interaction energy obtained from the 
Cq term can be decomposed into pair-wise contributions as 

A B C 

E d i S pfi = ^2 ( 6 ) 

a b Uab 

with 

a a a b (Mt) a (Mp b 
6 ' a6 a a (M?) b + a b (M?) a ^ 

In the expressions for these interatomic coefficients the polarizabilities and expec- 
tation values of the squares of the multipole moments of the systems are replaced 
by atomic polarizabilities and expectation values of the squares of the atomic multi- 
pole moments. In our previous worlP-2 we suggested an all-round Hirshfeld approach, 
where the Hirshfeld atomic multipole moments^Eni anc j Hirshfeld atomic intrinsic 
polarizabilities^ are employed. 



The Hirshfeld method allows to partition properties into atomic contributions by 
means of a weight function. The weight of each atom is determined by the density of 
the corresponding free spherical atom, normalized by the sum of all the free atomic 
densities of the atoms in the molecule 

freer \ 

UA (r) = W , (8) 

The elements of the intrinsic atomic polarizability tensor are then defined by^U 

<4 = j uw^(r)p^(r)dr, (9) 

where % and j represent the Cartesian directions x, y or z and pV' (r) denotes the 
first order density perturbed by an electric field applied in direction j. Since the 
polarizability is not a straightforwardly additive property, the total polarizability of 
the molecule cannot be reconstructed from the intrinsic polarizabilities alone, but 
a charge transfer term must be added. However, in the present work we will only 
consider the intrinsic polarizabilities. 

In this work the method is further extended by introducing the improved, itera- 
tive Hirshfeld approach (Hirshfeld-I), recently developed by Bultinck et alP^ to the 
calculation of dispersion coefficients. The Hirshfeld-I method differs from the clas- 
sic Hirshfeld method (Hirshfeld-C) by the definition of the atomic weight function. 
Whereas in Hirshfeld-C the weight function is predefined by the atomic densities of 
the free spherical atoms, in Hirshfeld-I the weight function is iterated until self consis- 
tency and therefore loses its somewhat arbitrary character. The Hirshfeld-C weight 



function is used as a first guess, and in each consecutive iteration the new weight 
function is constructed from the atomic densities obtained from the weight function 
of the previous iteration 

^sSw (10) 

The process is repeated untill the weight functions of two subsequent iterations are 
identical. 

The use of atomic intrinsic polarizabilites, which are obtained in the form of an 
atomic polarizability tensor, allows to introduce anisotropy into the model described 
above. 

Going back to the classical paper by Buckingham,®^ standard second order 
Rayleigh-Schrodinger perturbation theory for long range intermolecular forces was 
shown to yield the following general expression for the i?~ 6 contribution to the dis- 
persion energy for two molecules or two atoms a and b 

E abtdisp (R-°) = - ^f + b Uh) E E E E T^T^af (11) 

where aS a > and are the (dipole) polarizability tensors of the interacting systems 
and the elements of the T 2 tensor are defined as 

T 2 ,ij = ViVjR-, 1 (12) 

where R a b is the intermolecular distance and i and j stand for the Cartesian coordi- 
nates x,y,z. This expression was obtained using an Unsold/London type of approx- 
imatiorP^EH by simplifying the energy denominator in the second order perturbation 
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theory expression by introducing an average excitation energy (or ionization energy), 
for a and 6, namely U a and U b . Equation ffl~T]) can then be rewritten as 



E ab , disp (R- 6 ) = - 4{Ua + Ub) Tr (Ua } U a ) (13) 

<-> 

where the n and a tensors involve the elements and aiy. It is easily shown that 
in the case of isotropic tensors equation ([TBI reduces to 

E ab , dlsp (R-«) = -\ ±- -ML. a (a) a 9) Tr (fi 2 ) (14) 
4 K ab U a + tVb 

where a^ a > and cr ^ are now equal to the diagonal elements of a and a (or one 
third of their traces), respectively. As it is easily shown that Tr (n ) = 6, eq. (fl3|) 



finally reduces to 



whereby the standard London dispersion formula is recovered. Concentrating now on 
a pairwise atom-atom interaction scheme, where in (fI3"|) a and b refer to isolated atoms 
or atoms-in-molecules, and replacing in (II 3ft. in the spirit of Becke and Johnson's 
treatment!^ as also adopted in our previous work, the average excitation energies U a 
and Ub by expressions of the type 2 (Mf) /3a, where a is the isotropic polarizability 
of the atom in the molecule, we arrive at 

1 (Mf) (a) (Mf) (b) ^(a)^b 

G<Xb(Mf) {a) + a a {M{) {b) 
In the isotropic case the equation reduces to 

E (R-«\ - {Ml) ^ {Ml) ^ a a ± (17) 

dlSpM } * b (M?) {a) + a a (M?) (b) aaab Rl b [U) 
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which is the expression for the R~ 6 term used in our previous workP^ 

Equation ffT6|) was implemented in the atdisf^ program. Note that anistropy 
corrections to the R~ 8 and R~ w terms could be treated in a similar way necessitating, 
however, the evaluation of terms involving quadrupole and mixed dipole-quadrupole 
polarizabilities. As these contributions are expected to be smaller and lend themselves 
less easily to an interpolation into the framework of eq. ffT6l) . we only consider in 
this paper the expression of an anisotropy corrected C 6 term, optionally combined 
with isotropic C% and C\q terms. On the whole, anisotropy corrections to dispersion 
coefficients were relatively seldom studied in the literature.^ 

4 Computational Details 

The main goal of this work is to reproduce accurate interaction energies obtained from 
high level calculations by adding dispersion energy corrections to interaction energies 
obtained at the DFT level. For this purpose, a set of eighteen different complexes was 
examined. In order to ensure that the benchmark set is of good and consistent qual- 
ity, individual geometries of the different complexes have been optimized as follows. 
First the geometries of the monomers have been fully optimized at the CCSD(T)/aug- 
cc-pVTZ level with tight convergence criteria. Subsequently, the geometries of the 
complexes were optimized at the same level keeping the internal geometries of the 
monomers frozen. For most of the complexes this meant optimizing only one param- 
eter, namely the distance between the two monomers. For a few others the lateral 
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displacement of the two monomers has been taken into account. Symmetries of the 
complexes with the lowest reported energies were taken from the literature.^^^^ 
Once the equilibrium geometries were obtained, the interaction energies at CCSD(T) 
level and at DFT level using the B3LYP, PBE and TPSS functionals were calculated, 
taking into account the counterpoise BSSE corrections^ Since the geometry of the 
monomers was kept unchanged in the dimers, the interaction energies are given by 



W CCSD{T) _ rpCCSD(T) 



E 



CCSD(T) 



E 



CCSD(T) 
B 



rpDFT rpDFT f rpDFTl * f rpDFTl ' 
Winter = h 'AB ~ [ h A J ~ [ h B \ 



18) 



(19) 



The stars in eq. ffTSl) and fTTTJT) denote that the energies were obtained using all the 
basis functions of the dimer. 

To ensure high quality reference data the values for the interaction energy at the 
CCSD(T) level were extrapolated to the complete basis set limit (CBS) using the 
following focal point analysis. To obtain the BSSE-corrected interaction energy one 
needs to extrapolate the energies of the dimers and the complexes 



p CCSD(T)/CBS _ j^CCSD{T)/CBS 



E 



CCSD(T)/CBS 



E 



CCSD{T)/CBS 



B 



(20) 



where the total energy of each entity is a sum of the Hartree-Fock energy and the 



correlation energy 



^CCSD(T)/CBS _ pHF/CBS W CCSD{T)/CBS 



J tot 



(21) 



Since it has been shown by Sinnokrot and Sherrilp2] that the correlation energies at 
the MP2 and CCSD(T) levels converge very similarly with the size of the basis set, 
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the difference between the two remaining constant for the set of the aug-cc-pVXZ 
basis sets upon the increase of X, it is sufficient to extrapolate the correlation energy 
at the MP2 level to the basis set limit and add the correlation energy value between 
the two methods obtained for a smaller basis set: 

rpCCSD(T)/CBS _ rpMP2/CBS , r pCCSD(T)/aug-cc-pVTZ _ pMP2/aug-cc-pVTZi tcyy\ 
corr corr ' I corr corr J \ ) 

The MP2 correlation energy at the complete basis set limit was obtained using Hel- 
gaker's linear extrapolation formulaP^' 

-^3 rpMP2/aug-cc-pVXZ v 3 rpMP2/aug-cc-pVYZ 

T?MP2/CBS _ A &corr — Y theory 

corr 

where we chose X=5 and Y=4 for all the complexes. Finally, since the Hartree-Fock 
total energy converges towards the complete basis set limit fast and monotonously,^ 
the energies at the HF/aug-cc-pV6Z level were used to estimate Ef Q f S . 

In the following step the interaction energy at the DFT level, obtained with the 
B3LYP, PBE and TPSS functionals, was corrected for dispersion 

rpcorr.DFT _ rpDFT , rpDFT (r>A\ 
^inter ~ ^inter ^disp ' l z V 

where the dispersion interaction correction was calculated in four ways of increasing 
complexity. The first two approaches, which neglect anistropy, yield the following 
equations: 

• Isotropic C*6 term only 

A B s-iiso 

e dft _ rrV (2 *\ 

a b ah 
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Three isotropic terms 

rpDFT _ I ^6,ab w 8,afe ^lO.afr \ ( 0Pi \ 

^dispjull - / A nfi + + pin I \ £Ki ) 

^ab ^ab ^ab 



A B i r-iiso rnso r~f 



ISO 



a b 



In a subsequent step anisotropy is introduced. In order to reach an optimal quality 
cost ratio only the C$ term was corrected for anistropy as it may be expected that 
the main contribution to anisotropy will essentially be due to the C$ term. The 
counterpart of eq. f l25j) may be written as 

• Anisotropic C§ term only 

A B 

rpDFT _ _ \ r \ * rpaniso / r>— 6\ (97} 

^dispfi^ 130 ~ / j / j Cj disv,ab\ IX ) \ z ' J 



a b 



• A compilation of the anisotropic C 6 term and isotropic C 8 and C w terms finally 
gives 

A B / /~iiso riiso \ 

rpDFT _ _ I paniso I r>-6\ , 8,afc °10,ab \ 

11 'disp, mixed ~ / 4 / 4 I - C/ disp,afel- rl I ' p8 ' DlO I ' V / 

a 6 \ ^ ^afc ) 

The coefficients used in eq. (I25p to (|28p were obtained by the Hirshfeld-I method. Note 
that no damping function has been used in these equation, since we are interested 
here in dispersion energies of complexes at their equilibrium geometry, where the 
interatomic distances are relatively large. 

The geometries and polarizabilities of the different molecules were calculated us- 
ing the GAUSSIAN0 3 65 program. The values of atomic intrinsic polarizabilities and 
expectation values of the squared atomic multipole moments were obtained using the 
g TOCK |39] p ro g ram Finally the different dispersion coefficients and the corresponding 
dispersion energy corrections were calculated using the atdisp^ program. 
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5 Results and Discussion 

Table 1 lists the interaction energies obtained with the CCSD(T)/CBS, B3LYP/aug- 
cc-pVTZ, PBE/aug-cc-pVTZ and TPSS/aug-cc-pVTZ methods. The high level in- 
teraction energies vary from -0.02 to -1.76 kcal/mol, whereas the interaction energies 
obtained at the DFT level seem to be dependent on the choice of functional. For 
the B3LYP functional all interaction energy values are repulsive, as expected at the 
DFT level where non-local dispersion interactions are not included in the exchange- 
correlation functional. However, PBE produces negative interaction energies for all 
of the examined complexes while the values obtained with the TPSS functional are 
partially negative. Since the attraction for most of the examined complexes can be 
attributed only to dispersion energy, the negative interaction values at the DFT level 
are not physically justified, B3LYP thus being in this aspect the most reliable of the 
three functionals examined. 

Tables 2 to 4 give the dispersion corrected post-DFT interaction energies obtained 
by the C^-isotropic (eq. [25j) , full isotropic (eq. [26j) , Cq anisotropic (eq. [2TJ) and mixed 
(eq. I2"8j) models for the B3LYP, PBE and TPSS functionals, respectively. In contrast 
to our previous work on dispersion energies,^ where the classic version of the Hirshfeld 
method was utilized, the dispersion corrections are obtained here with the Hirshfeld-I 
method. The difference between the two versions of the method is only of importance 
for the larger complexes, where the monomers contain more than one atom, which are 
unidentical. In that case the atomic weight functions which determine the distribution 
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of the electronic density among the atoms are different, and as a result also the 
charges and other properties of the atoms are different. It has been shown in a 
previous study by some of the authors^ that atomic polarizabilities obtained with 
the Hirshfeld-I method are of better quality, therefore also leading to more reliable 
dispersion coefficients. The two last lines in Tables 2 to 4 also show the correlation 
coefficient between the post-DFT interaction energies and the high level interaction 
energies listed in Table 1 and the standard error of the linear regression. 

From Table 2 it appears that for the B3LYP functional, the addition of a disper- 
sion interaction correction based on the Cq dispersion coefficient reduces most of the 
interaction energies to negative values, although for five complexes the interaction 
energies remain positive. Taking into account the remaining Cg and C±o coefficients 
further reduces all post-B3LYP interaction energies to values similar to the high-level 
interaction energies. Figure 1 depicts the correlation between the post-B3LYP in- 
teraction energy values, obtained using the isotropic and anisotropic C 6 dispersion 
coefficients, and the high level values. Although the correlation is poor, being only 
R = 0.7064 and R = 0.8098 for the isotropic and anisotropic models, respectively, 
the effect of the addition of anisotropy can be seen very clearly in this Figure. Obvi- 
ously, the anisotropic contributions are only of significance for the complexes where 
the monomers are not spherically symmetric, so the values for the smaller complexes 
remain the same. While for the isotropic model the values for the CO2-CO2 and C2H4- 
C2H4 complexes are situated far from the trend line, introduction of the anisotropy to 
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the C 6 coefficient lowers the value for the C 2 H 4 -C 2 H4, although the C0 2 -C0 2 complex 
remains an outlier. One can conclude that the large standard error (listed in Table 2) 
and the very low value of the slope (shown in Figure 1), indicate that the CVbased 
models are insufficient by far for reproducing accurate post-DFT interaction energies 
and that the higher coefficients are indispensable. 

Figure 2 depicts the correlation between the high level interaction energy values 
and the post-B3LYP interaciton energy values for the full isotropic and mixed models. 
As was mentioned in the method section, the derivation of anisotropic dispersion co- 
efficients higher than C 6 becomes complicated as extra terms appear in the equations 
and the model loses its simplicity, which is one of its major advantages. Therefore 
a mixed model is examined here, which on the one hand combines the improvement 
achieved in the C 6 coefficient by introducing anisotropy and on the other hand uses 
isotropic Cg and Cio coefficients, which are vital for reproduction of accurate interac- 
tion energies. The mixed model performs better than the full isotropic method, both 
of them performing significantly better than the CVbased models, having a correla- 
tion coefficient of 0.94 and 0.95 for the former and the latter, respectively. Also the 
value for the C0 2 -C0 2 complex improves considerably by the addition of the higher 
dispersion coefficients. The difference between the two models is the largest for the 
four complexes with the highest interaction energy values, namely C0 2 -C0 2 , OCS- 
OCS, C 2 H 2 -C 2 H 2 and C 2 H 4 -C 2 H 4 . As can be seen in Figure 2, when only isotropic 
contributions are taken into account, the order of interaction energies for those three 
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complexes is not correctly reproduced: while CCSD(T) interaction energies follow the 
order (in absolute values) OCS-OCS > C2H4-C2H4 > C 2 H 2 -C 2 H 2 > C0 2 -C0 2 , the in- 
teraction energies obtained with the full isotropic model follow the order C 2 H 4 -C 2 H 4 
> OCS-OCS > C 2 H 2 -C 2 H 2 > C0 2 -C0 2 . However, once anisotropy is introduced 
in the mixed model, the correct order is restored. The reason for this difference is 
the lower symmetry of the complexes, where the monomers are laterally shifted with 
respect to each other and anisotropy becomes more important. Therefore one can 
expect the linear regression parameters to improve further for the mixed model when 
more complexes of lower symmetry are taken into account. 

Although the correlation achieved for the full isotropic and mixed models is sat- 
isfying, the linear regression parameters are not optimal yet. With a slope of 0.82 in 
the full isotropic model and 0.83 in the mixed model, both models underestimate the 
interaction energy by an approximate 20%. The PBE functional seems to suffer from 
an opposite problem, as can be seen from the values listed in Table 3. The interaction 
energies obtained by the full and mixed models are overestimated by no less than 40%, 
as can also be seen from the linear regression coefficients in Figure 3. The source of 
the overestimation of the values must be sought in the pure DFT interaction energies 
obtained with this functional. As can be seen in Table 1, the PBE functional pro- 
duces all negative interaction energies even though dispersion is not included in the 
functional. Since the dispersion energy correction is added to the original pure DFT 
interaction value, the spurious potential well produced by the PBE functional for the 
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examined complexes causes a serious overestimation of the interaction energies. On 
the other hand, the correlation coefficients for the full isotropic and mixed models 
are very high, being more than 0.99, while the standard error is reduced almost half 
in size in comparison to the B3LYP functional. One can therefore conclude that our 
method for the calculation of dispersion energy corrections performs surprisingly well 
for the PBE functional, but for accurate interaction energies the values must be scaled 
by a factor of 0.71. It must also be mentioned that the effect of anisotropy on the 
interaction energy values is analogous for the PBE functional: the replacement of the 
isotropic C 6 coefficient by an anisotropic one in the mixed model restores the correct 
order in interaction energy values for the four largest complexes. The addition of the 
higher order coefficients C 8 and C w is here also of importance. Although the corre- 
lation coefficients are quite high for the C 6 -based models, the values are less reliable. 
For example, the interaction value for the C2H2-C2H2 is an evident outlier in Table 3 
for those two models. 

The post-TPSS interaction energies are listed in Table 4 for the four different 
models and depicted in Figure 4 for the full isotrpic and mixed models. This functional 
seems to perform very well, the correlation for the full isotropic and mixed models 
being 0.98 and the standard error on the linear regression being almost as low as 
for the PBE functional. The main strength of this functional is the perfect slope 
of 1.00 for the mixed model, enabling us to reproduce accurate interaction energies 
without the need for up- or down-scaling, as is the case for the B3LYP and PBE 
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functionals. Amongst the larger complexes only the interaction energy value of the 
ethene dimer appears to be problematic, being overestimated by 0.4 kcal/mol. As 
was also the case with the PBE functional, the source for this overestimation may lay 
in the negative pure DFT interaction energy value of this complex (-0.548 kcal/mol). 
It seems that also here our method for obtaining dispersion energies is performing 
very well, but care must be taken if the potential energy surface produced by the 
functional is incorrect. 



6 Concluding Remarks 

In our previous workpH it was shown that the combination of Becke and Johnson's 
exchange dipole moment approach and our Hirshf eld-type partitioning scheme for 
molecular polarizabilities lead to a simple, accurate and inexpensive approach for 
evaluation of dispersion energies of dimers. This approach was further developed to 
reproduce the chemically more interesting interaction energies at the CCSD(T) level 
with simple dispersion energy corrected DFT interaction energies. Three function- 
als different by nature were examined here to test the robustness of our proposed 
method, namely B3LYP, PBE and TPSS. Two major conclusions can be drawn from 
the present results. First of all, the inclusion of anisotropy in the evaluation of the 
Cq coefficient leads to an improvement of the corresponding dispersion energies. This 
could be seen especially for the dispersion energy values of the four complexes with 
lower symmetry, namely C0 2 -C0 2 , OCS-OCS, C 2 H 2 -C 2 H2 and C 2 H 4 -C 2 H 4 , where 
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the mixed model produced significant improvements. Secondly, our method performs 
well for different functionals regardless of their nature and might therefore be uni- 
versally applicable to different DFT functionals. However, since the final result for 
the interaction energy is not only dependent on our dispersion energy correction but 
also on the pure DFT interaction energy, care must be taken in the choice of the 
functional. From the three functionals examined here TPSS seems to perform the 
best, producing accurate interaction energies for most of the complexes. However, for 
complexes where the TPSS functional produces questionable pure DFT interaction 
energies, such as for example was the case for the C2H2-C2H2 complex, the B3LYP 
functional is a better choice. 
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8 Tables 



Complex 


CCSD(T) 


B3LYP 


PBE 


TPSS 


He-He 


-0.020 


0.041 


-0. 


.056 


-0.043 


He-Ne 


-0.037 


0.040 


-0. 


083 


-0.058 


He-Ar 


-0.056 


0.072 


-0. 


.078 


-0.047 


Ne-Ne 


-0.071 


0.045 


-0. 


108 


-0.066 


Ne-Ar 


-0.119 


0.021 


-0. 


123 


-0.056 


Ar-Ar 


-0.272 


0.177 


-0. 


120 


0.017 


L-He-N2 


-0.043 


0.088 


-0. 


.073 


-0.033 


T-He-N2 


-0.061 


0.109 


-0. 


.075 


-0.027 


He-FCl 


-0.096 


0.072 


-0. 


116 


-0.046 


FCl-He 


-0.126 


0.048 


-0. 


.220 


-0.079 


Ne-CH 4 


-0.175 


0.066 


-0. 


.123 


0.022 


CH4-C2H4 


-0.449 


0.428 


-0. 


101 


0.158 


CH4-CH4 


-0.537 


0.490 


-0. 


.025 


0.273 


SiH4-CH4 


-0.824 


0.549 


-0. 


.130 


0.312 


C2H2-C2H2 


-1.403 


0.153 


-0. 


.928 


-0.548 


CO2-CO2 


-1.476 


0.757 


-0. 


.388 


0.133 


OCS-OCS 


-1.761 


0.761 


-0. 


.083 


0.572 


C2H4-C2H4 


-1.493 


0.544 


-0. 


.284 


0.431 



Table 1: Interaction energies calculated with CCSD(T)/CBS, B3LYP/aug-cc-pVTZ, 
PBE/aug-cc-pVTZ and TPSS/aug-cc-pVTZ methods. All values are in kcal/mol. 



29 



Complex 


r^tso 
U 6 


r^aniso 


full 


mixed 


He-He 


0.01 


0.01 


-0.01 


-0.01 


He-Ne 


-0.01 


-0.01 


-0.04 


-0.04 


He-Ar 


0.01 


0.01 


-0.04 


-0.04 


Ne-Ne 


-0.03 


-0.03 


-0.08 


-0.08 


Ne-Ar 


-0.09 


-0.09 


-0.19 


-0.19 


Ar-Ar 


-0.08 


-0.08 


-0.34 


-0.34 


L-He-N2 


0.02 


0.03 


-0.02 


-0.01 


T-He-N2 


0.02 


0.01 


-0.05 


-0.05 


He-FCl 


-0.04 


-0.04 


-0.13 


-0.13 


FCl-He 


-0.04 


-0.04 


-0.09 


-0.09 


Ne-CH 4 


-0.15 


-0.15 


-0.28 


-0.28 


CH4-C2H4 


0.06 


0.03 


-0.11 


-0.14 


CH4-CH4 


-0.06 


-0.06 


-0.30 


-0.30 


SiH 4 -CH 4 


-0.21 


-0.19 


-0.62 


-0.61 


C2H2-C2H2 


-0.50 


-0.48 


-1.16 


-1.15 


CO2-CO2 


0.03 


-0.04 


-0.75 


-0.81 


OCS-OCS 


-0.33 


-0.47 


-1.56 


-1.71 


C2H4-C2H 4 


-0.69 


-0.55 


-1.64 


-1.50 


R 


0.7064 


0.8098 


0.9412 


0.9544 


a 


0.44 


0.37 


0.18 


0.16 



Table 2: The four different types of post-DFT interaction energies calculated with 
the Hirshfeld-I method and the B3LYP functional. All values are in kcal/mol. 



30 



Complex 


r^tso 
^6 


r^aniso 


full 


mixed 


He-He 


-0.09 


-0.09 


-0.11 


-0.11 


He-Ne 


-0.13 


-0.13 


-0.17 


-0.17 


He-Ar 


-0.15 


-0.15 


-0.20 


-0.20 


Ne-Ne 


-0.19 


-0.19 


-0.25 


-0.25 


Ne-Ar 


-0.24 


-0.24 


-0.35 


-0.35 


Ar-Ar 


-0.38 


-0.38 


-0.64 


-0.64 


L-He-N2 


-0.14 


-0.13 


-0.18 


-0.18 


T-He-N2 


-0.17 


-0.17 


-0.24 


-0.24 


He-FCl 


-0.23 


-0.24 


-0.33 


-0.33 


FCl-He 


-0.31 


-0.31 


-0.37 


-0.37 


Ne-CH 4 


-0.35 


-0.35 


-0.49 


-0.49 


CH4-C2H4 


-0.48 


-0.54 


-0.66 


-0.71 


CH4-CH4 


-0.59 


-0.57 


-0.85 


-0.82 


SiH 4 -CH 4 


-0.92 


-0.90 


-1.35 


-1.34 


C2H2-C2H2 


-1.58 


-1.57 


-2.25 


-2.24 


CO2-CO2 


-1.15 


-1.21 


-1.96 


-2.02 


OCS-OCS 


-1.19 


-1.33 


-2.44 


-2.59 


C2H4-C2H 4 


-1.54 


-1.40 


-2.50 


-2.36 


R 


0.9595 


0.9777 


0.9902 


0.9953 


a 


0.18 


0.13 


0.12 


0.09 



Table 3: The four different types of post-DFT interaction energies calculated with 
the Hirshfeld-I method and the PBE functional. All values are in kcal/mol. 
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Complex Q so C% niso full mixed 



He-He 


-0. 


.07 


-0 


,07 


-0 


.09 


-0 


,09 


He-Ne 


-0 


.11 


-0. 


,11 


-0, 


.14 


-0 


.14 


He-Ar 


-0. 


.11 


-0 


.11 


-0. 


.16 


-0, 


.16 


Ne-Ne 


-0 


.14 


-0 


.14 


-0. 


.20 


-0. 


.20 


Ne-Ar 


-0. 


.17 


-0 


.17 


-0 


.28 


-0 


.28 


Ar-Ar 


-0 


.24 


-0. 


.24 


-0 


.49 


-0 


.49 


L-He-N2 


-0 


.10 


-0. 


.09 


-0 


.14 


-0 


.13 


T-He-N2 


-0 


.12 


-0 


.12 


-0 


.18 


-0 


.18 


He-FCl 


-0. 


.16 


-0. 


.16 


-0 


.25 


-0. 


,25 


FCl-He 


-0. 


.16 


-0. 


.16 


-0 


.22 


-0 


,22 


Ne-CH 4 


-0 


.20 


-0. 


.19 


-0. 


.33 


-0. 


.33 


CH4-C2H4 


-0. 


.21 


-0. 


,27 


-0. 


.37 


-0, 


,43 


CH4-CH4 


-0 


.27 


-0. 


.25 


-0. 


.51 


-0. 


,48 


SiH 4 -CH 4 


-0 


.43 


-0. 


,42 


-0 


.83 


-0. 


,82 


C2H2-C2H2 


-1 


.19 


-1. 


.18 


-1, 


.83 


-1. 


.82 


CO2-CO2 


-0 


.61 


-0 


.67 


-1, 


.38 


-1, 


.45 


OCS-OCS 


-0 


.50 


-0 


,65 


-1 


.70 


-1, 


.85 


C2H4-C2H 4 


-0 


.79 


-0 


,65 


-1 


.72 


-1, 


.58 



R 0.8572 0.8934 0.9791 0.9853 

a 0.32 0.28 0.13 0.11 



Table 4: The four different types of post-DFT interaction energies calculated with 
the Hirshfeld-I method and the TPSS functional. All values are in kcal/mol. 
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9 Figures 
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Figure 1: Correlation between high level interaction energies and post-B3LYP inter- 
action energies obtained with the isotropic and anisotropic C 6 models. All values are 
in kcal/mol. 
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Figure 2: Correlation between high level interaction energies and post-B3LYP inter- 
action energies obtained with the full isotropic and mixed models. All values are in 
kcal/mol. 



34 



CCSD(T) interaction energy 

-2.3 -1.8 -1.3 -0.8 -0.3 0.2 

i 1 1 1 1 1 0.2 




♦ Full isotropic □ Mixed 



Figure 3: Correlation between high level interaction energies and post-PBE interac- 
tion energies obtained with the full isotropic and mixed models. All values are in 
kcal/mol. 
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Figure 4: Correlation between high level interaction energies and post-TPSS inter- 
action energies obtained with the full isotropic and mixed models. All values are in 
kcal/mol. 
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